(*
 *******************************************************
	Package: Phase2D
	Author:  Janos Karsai
	Address: Dept. Medical Informatics,
		 University of Szeged
	         H-6720, Szeged, Hungary
	E-mail:  karsai@dmi.u-szeged.hu
 *******************************************************
*)

(*
Package for visualization of solutions of families of planar parametric curves
version: Mathematica 3.0
*)


(*
*******************************************************
Developed for the book -CDRom: 

	Models of Impulsive systems

	Typotex, Budapest, 2002

* Development is supported by 
* Foundation for the Hungarian Higher Education Research, Grant No. 776/94
* Foundation for the Development of Higher Education Grant No. FEFA-IV 1608/4.
* Hungarian Science Research Foundation OTKA T-034275 
* Hungarian Ministry of Education: Project o. 179/1996
* Pro Renovanda Cultura Hungarie: "Science in Education"
* Hungarian Ministry of Education: Higher Education Textbook Grant, 2001
*
*******************************************************
*)


(* COPYRIGHT_BEGIN
*  Copyright (c) 2002 János Karsai,
*  Department of Medical Informatics, University of Szeged
*
*    
*  Permission to use, copy, modify, distribute, and sell this software and its documentation
*  for any purpose is hereby granted without fee, provided that (i) the above copyright 
*  notices and this permission notice appear in all copies of the software and related
*  documentation, and (ii) the name J. Karsai may not be used in 
*  any advertising or publicity relating to the software without the specific, prior written 
*  permission of J. Karsai.
*    
*  THE SOFTWARE IS PROVIDED "AS-IS" AND WITHOUT WARRANTY OF ANY
*  KIND, EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT 
*  LIMITATION, ANY WARRANTY OF MERCHANTABILITY OR FITNESS FOR
*  A PARTICULAR PURPOSE.
*    
*  IN NO EVENT SHALL J. KARSAI BE LIABLE FOR ANY SPECIAL,
*  INCIDENTAL, INDIRECT OR CONSEQUENTIAL DAMAGES OF ANY KIND,
*  OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA
*  OR PROFITS, WHETHER OR NOT ADVISED OF THE POSSIBILITY OF 
*  DAMAGE, AND ON ANY THEORY OF LIABILITY, ARISING OUT OF OR IN 
*  CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
* COPYRIGHT_END 
*)



(*
The functions
The functions of this package provide several 2D and 3D  plots of families of 2D trajectories, 
such as  trajectories of the solutions of the 2D system of differential equations:
		x'=xdot
		y'=ydot
with a list initial conditions.
*)
(*
The parameters
t0		:initial time
T		:final time
Sol		:{x[t],y[t]} a curve in the 2D phaseplain
Traj[t]		:List of trajectories of form {{x[t],y[t]},...} as functions
Traj3D[t]	:List of lists {t,x[t],y[t]}
*)

UnsortedComplement[x_List, y__List] := Replace[x, (# :> Sequence[]) & /@ Join[y], 1];

Clear[ParametricPlotColor]
Options[ParametricPlotColor]=Options[ParametricPlot];
Attributes[ParametricPlotColor] = {HoldAll};
automatic[x_, value_] := If[x === Automatic, value, x]; 
ParametricPlotColor[xy_List,{t_,t0_,t1_},opt___]:=
	Module[{color,xy2,xy2a,ddt,options,pltsty,pltopt,nonopt},
	
	options=List[opt,Sequence@@Options[ParametricPlotColor]];
	
	ddt = automatic[(t1 - t0)/PlotPoints/.options, (t1 - t0)/30]; 
	pltsty=automatic[PlotStyle/.options, {}]; 
		nonopt={PlotPoints,PlotDivision,Compiled,MaxBend,PlotStyle};
	
(*
pltopt=UnsortedComplement[ 
        options /.(
            nonopt /.a_Symbol->((a->b_)->a)/.((List->b_)->List) ->List),nonopt];	
	pltopt[[0]]=Sequence;
*)
	pltopt=UnsortedComplement[
  		Map[Which[MemberQ[nonopt, #[[1]]], KJ, True, #]&, options], {KJ}];

Return[
		Which[Length[Dimensions[xy]]==1 && Length[xy]==2,
			{ xy2=xy /.t->u; xy2a=xy2/.(u->u+ddt);
			Show[Graphics[
				Table[{pltsty,
					Line[{xy2,xy2a}]},
					{u,t0,t1-ddt,ddt}]],
      				Sequence@@pltopt]},
		      Length[Dimensions[xy]]==1 && Length[xy]==3,
			{color=xy[[3]]/. t:>u; xy2=Take[xy,2]/.t:>u; xy2a=xy2/.(u->u+ddt);
			Show[Graphics[
				Table[{pltsty,color,
             				Line[{xy2,xy2a}]},
					{u,t0,t1-ddt,ddt}]],
				Sequence@@pltopt]},
		      True,
			Show[
				Map[Graphics[
					If[Length[#]==2,
					Table[ Join[ pltsty,{Line[{#/.t->u,(#/.t->u+ddt)}]}],
						{ u,t0,t1-ddt,ddt}],
					Table[Join[ pltsty,{#[[3]] /.t->u,
                     				Line[{Drop[#,-1]/.t->u,(Drop[#,-1]/.t->u+ddt)}]}],
						{ u,t0,t1-ddt,ddt}]]
				    ]&,xy],
				Sequence@@pltopt] 
		]
	];
	];


(* 2D plots *)

(* Plot the trajectories *)
PhasePlot[Traj_List,{t_,t0_,t1_},params___]:=
	ParametricPlotColor[Evaluate[Map[Append[#,Hue[(t-t0)/(t1-t0)]]&,Traj]],
		{t,t0,t1},params, AspectRatio->Automatic,PlotStyle -> {Thickness[0.015]}];

(* Plot the trajectories *)
PhasePlotBW[Traj_List,{t_,t0_,t1_},params___]:=
	ParametricPlot[Evaluate[Map[(#)&,Traj]],
		{t,t0,t1},params, AspectRatio->Automatic,PlotStyle -> {Thickness[0.015]}];

ListPhasePlot[sol_,lineopt_List, opt___] := 
  Show[Map[ListPlot[Coordinate[#, {2, 3}], PlotStyle->Prepend[lineopt,Thickness[0.015]],
 	  DisplayFunction -> Identity, 
          PlotJoined -> True,
          AspectRatio -> Automatic] &, sol], opt, 
    DisplayFunction -> $DisplayFunction]

ListPhasePlot[sol_, opt___] := 
  Show[Map[ListPlot[Coordinate[#, {2, 3}], DisplayFunction -> Identity, 
          PlotJoined -> True, PlotStyle -> {Thickness[0.015]}, 
          AspectRatio -> Automatic] &, sol], opt, 
    DisplayFunction -> $DisplayFunction]


(* The phaseMap at t1 *)
PhaseMap[Traj_,{t_,t1_},opt___]:=Show[Graphics[{Line[Traj /.t->t1]}],
		opt,
		AspectRatio->Automatic, 
		PlotLabel->StringJoin["t=",ToString[N[t1]]]];

PhaseMap[Traj_,{t_,t1_},lineopt_List,opt___]:=
		Show[Graphics[Flatten[{lineopt,Line[Traj /.t->t1]}]],
		opt,
		AspectRatio->Automatic, 
		PlotLabel->StringJoin["t=",ToString[N[t1]]]];

PhaseMap[Traj_,{t_,t0_,t1_,dt_},opt___]:=		
		Table[
		Show[
		Graphics[Flatten[{Thickness[0.02],Hue[(s-t0)/(t1-t0)],Line[Traj /.t->s]}]],
		opt,AspectRatio->Automatic, 
		PlotLabel->StringJoin["t=",ToString[N[s]]]],{s,t0,t1,dt}];

PhaseMap[Traj_,{t_,t0_,t1_,dt_},lineopt_List,opt___]:=		
		Table[
		Show[
		Graphics[Flatten[{Thickness[0.02],Hue[(s-t0)/(t1-t0)],lineopt,Line[Traj /.t->s]}]],
		opt,AspectRatio->Automatic, 
		PlotLabel->StringJoin["t=",ToString[N[s]]]],{s,t0,t1,dt}];

(* The phaseMap at t1 from the left and right *)



PhaseMapImp[Traj_,{t_,t1_?NumberQ},lineopt_List,opt___]:=
		Show[Graphics[
		Flatten[{Dashing[{0.1,0.03}],lineopt,Line[Traj /.t->(t1-0.001)]}]],
		Graphics[Flatten[{lineopt,Line[Traj /.t->t1]}]],
		opt,
		AspectRatio->Automatic, 
		PlotLabel->StringJoin["t=",ToString[N[t1]]]];

PhaseMapImp[Traj_,{t_,t1_?NumberQ},opt___]:=Show[
		Graphics[{Dashing[{0.1,0.03}],Line[Traj /.t->(t1-0.001)]}],
		Graphics[{Line[Traj /.t->t1]}],
		opt,
		AspectRatio->Automatic, 
		PlotLabel->StringJoin["t=",ToString[N[t1]]]];

PhaseMapImp[Traj_,{t_,tn_List},lineopt_List,opt___]:=		
		Map[
		Show[
		Graphics[
		Flatten[{Hue[(#-First[tn])/(Last[tn]-First[tn])],Dashing[{0.03,0.03}],lineopt,
		Line[Traj /.t->(#-0.001)]}]
		],
		Graphics[Flatten[{Hue[(#-First[tn])/(Last[tn]-First[tn])],lineopt,Line[Traj /.t->#]}]
		],
		opt,AspectRatio->Automatic, 
		PlotLabel->StringJoin["t=",ToString[N[#]]]]&,tn];

PhaseMapImp[Traj_,{t_,tn_List},opt___]:=		
		Map[
		Show[Graphics[
		Flatten[{Hue[(#-First[tn])/(Last[tn]-First[tn])],Dashing[{0.03,0.03}],
			Line[Traj /.t->(#-0.001)]}]],
		Graphics[Flatten[{Hue[(#-First[tn])/(Last[tn]-First[tn])],Line[Traj /.t->#]}]],
		opt,AspectRatio->Automatic, 
		PlotLabel->StringJoin["t=",ToString[N[#]]]]&,tn];

(* The ListPlot of the phaseMap at t1 *)

ListPhaseMap[Traj_,{t_,t1_},opt___]:=ListPlot[Traj /.t->t1,
		opt,
		AspectRatio->Automatic, 
		PlotLabel->StringJoin["t=",ToString[N[t1]]]];

ListPhaseMap[Traj_,{t_,t0_,t1_,dt_},opt___]:=
		Table[ListPlot[Traj /.t->s,opt,
		AspectRatio->Automatic, PlotStyle -> {PointSize[0.04], Hue[s/(t1 - t0)]},
		PlotLabel->StringJoin["t=",ToString[N[t1]]]], {s, t0, t1, dt}];


(* Plot of the i-th coord. of the solutions *)

SolCoordPlot[Traj_,{t_,t0_,T_},Icoord_,opt___]:=
	(Plot[#,{t,t0,T},opt 
	 		 ] & )[Transpose[Traj][[Icoord]]]




(* Plots every soln. in Traj[t] in graph. array *)

Clear[SolPlot]
SolPlot[Traj_, {t_, t0_, T_}, Label_, opt_List] :=
  Map[Show[GraphicsArray[Partition[
            Table[Plot[Evaluate[#[[i]]], {t, t0, T}, PlotLabel -> Label[[i]]
                , DisplayFunction -> Identity, 
                Evaluate[Sequence @@ (opt[[i]])]], {i, 1, Length[#]}], 
            Length[#]]], DisplayFunction -> $DisplayFunction] &, Traj]

SolPlot[Traj_, {t_, t0_, T_}, Label_, opt___] :=
  Map[Show[GraphicsArray[Partition[
            Table[Plot[Evaluate[#[[i]]], {t, t0, T}, PlotLabel -> Label[[i]]
                , DisplayFunction -> Identity, opt], {i, 1, Length[#]}], 
            Length[#]]], DisplayFunction -> $DisplayFunction] &, Traj]


SolPlot[Traj_, {t_, t0_, T_}, opt___] :=
  Map[Show[GraphicsArray[Partition[
            Table[Plot[Evaluate[#[[i]]], {t, t0, T},
                DisplayFunction -> Identity, opt], {i, 1, Length[#]}], 
            Length[#]]], DisplayFunction -> $DisplayFunction] &, Traj]


Clear[ListSolPlot]

ListSolPlot[sol_, Label_, opt_List] := 
  Map[Show[GraphicsArray[
          Partition[
            Table[ListPlot[Coordinate[#, {1, i}], 
                Evaluate[Sequence @@ (opt[[i - 1]])], PlotJoined -> True, 
                PlotLabel -> Label[[i - 1]], DisplayFunction -> Identity], 
		{i, 2, Length[#[[1]]]}], Length[#[[1]]] - 1]], 
        DisplayFunction -> $DisplayFunction] &, sol]

ListSolPlot[sol_, Label_, opt_] := 
  Map[Show[GraphicsArray[
          Partition[
            Table[ListPlot[Coordinate[#, {1, i}], 
                opt, PlotJoined -> True, 
                PlotLabel -> Label[[i - 1]], DisplayFunction -> Identity], 
		{i, 2, Length[#[[1]]]}], Length[#[[1]]] - 1]], 
        DisplayFunction -> $DisplayFunction] &, sol]

ListSolPlot[sol_, opt_] := 
  Map[Show[GraphicsArray[
          Partition[
            Table[ListPlot[Coordinate[#, {1, i}], 
                opt, PlotJoined -> True, 
                DisplayFunction -> Identity], 
		{i, 2, Length[#[[1]]]}], Length[#[[1]]] - 1]], 
        DisplayFunction -> $DisplayFunction] &, sol]



(* ----------------------------------------------------- *)
(* 3D Plots *)
(* ----------------------------------------------------- *)

(* The solutions in the 3D space *)

Traj3D[Traj_,t_]:= Map[Join[{t},#]&,Traj];

(* Plot of the soln.-s in (t,x,y) system *)

SolPlot3D[Traj_,{t_,t0_, T_},opt___] := 
 ParametricPlot3D[Evaluate[Map[Append[#,Hue[(t-t0)/(T-t0)]]&,Traj3D[Traj,t]]],
		 {t, t0, T}, opt,
		AspectRatio -> Automatic, 
      		PlotRange -> All]

SolPlot3DBW[Traj_,{t_,t0_, T_},opt___] := 
 ParametricPlot3D[Evaluate[Map[(#)&,Traj3D[Traj,t]]],
		 {t, t0, T}, opt,
		AspectRatio -> Automatic, 
      		PlotRange -> All]

ListSolPlot3D[sol_, lineparm_, opt___] :=
  Show[Graphics3D[Join[{Thickness[0.01], RGBColor[0, 0, 0.8]},lineparm, Map[Line[#] &, sol]]], opt]


PhaseVol[Traj_,{t_,t1_},opt___]:=
	Show[	Graphics3D[Line[Traj3D[Traj,t]/.t->t1]],opt];

PhaseVol[Traj_,{t_,t1_},lineopt_List,opt___]:=
	Show[	Graphics3D[Flatten[{lineopt,Line[Traj3D[Traj,t]/.t->t1]}]],opt];


PhaseVol[Traj_,{t_,t0_,t1_,dt_},lineopt_List,opt___]:=				
Module[{showsol3d,s},
		showsol3d = SolPlot3D[Traj, {t, t0, t1},DisplayFunction->Identity];
		Table[
		Show[showsol3d,
		Graphics3D[
		Flatten[{Hue[(s-t0)/(t1-t0)],lineopt,
		Line[Traj3D[Traj,t] /.t->s]}]
		],
		opt,AspectRatio->Automatic, DisplayFunction -> $DisplayFunction,
		PlotLabel->StringJoin["t=",ToString[N[s]]]],{s,t0,t1,dt}]
	];

PhaseVolBW[Traj_,{t_,t0_,t1_,dt_},lineopt_List,opt___]:=				
Module[{showsol3d,s},
		showsol3d = SolPlot3DBW[Traj, {t, t0, t1},DisplayFunction->Identity];
		Table[
		Show[showsol3d,
		Graphics3D[
		Flatten[{Hue[(s-t0)/(t1-t0)],lineopt,
		Line[Traj3D[Traj,t] /.t->s]}]
		],
		opt,AspectRatio->Automatic, DisplayFunction -> $DisplayFunction,
		PlotLabel->StringJoin["t=",ToString[N[s]]]],{s,t0,t1,dt}]
	];
(* ----------------------- *)

PhaseVolImp[Traj_,{t_,t1_?NumberQ},lineopt_List,opt___]:=
	Show[	
	Graphics3D[Flatten[{lineopt,Dashing[{0.1,0.03}],Line[Traj3D[Traj,t]/. t->(t1-0.001)]}]],
	Graphics3D[Flatten[{lineopt,Line[Traj3D[Traj,t]/.t->t1]}]]
	,opt]; 


PhaseVolImp[Traj_,{t_,t1_?NumberQ},opt___]:=
	Show[	
		Graphics3D[{Dashing[{0.1,0.03}],Line[Traj3D[Traj,t]/.t->(t1-0.001)]}],	
		Graphics3D[Line[Traj3D[Traj,t]/.t->t1]]
	,opt];


PhaseVolImp[Traj_,{t_,tn_List},lineopt_List,opt___]:=				
Module[{showsol3d},
		showsol3d = SolPlot3D[Traj, {t, First[tn], Last[tn]},DisplayFunction->Identity];
		Map[
		Show[showsol3d,
		Graphics3D[
		Flatten[{Hue[(#-First[tn])/(Last[tn]-First[tn])],Dashing[{0.03,0.03}],lineopt,
		Line[Traj3D[Traj,t] /.t->(#-0.001)]}]
		],
		Graphics3D[Flatten[{Hue[(#-First[tn])/(Last[tn]-First[tn])],lineopt,
		Line[Traj3D[Traj,t]/.t->#]}]
		],
		opt,AspectRatio->Automatic, DisplayFunction -> $DisplayFunction,
		PlotLabel->StringJoin["t=",ToString[N[#]]]]&,tn]
	];

PhaseVolImpBW[Traj_,{t_,tn_List},lineopt_List,opt___]:=				
Module[{showsol3d},
		showsol3d = SolPlot3DBW[Traj, {t, First[tn], Last[tn]},DisplayFunction->Identity];
		Map[
		Show[showsol3d,
		Graphics3D[
		Flatten[{Hue[(#-First[tn])/(Last[tn]-First[tn])],Dashing[{0.03,0.03}],lineopt,
		Line[Traj3D[Traj,t] /.t->(#-0.001)]}]
		],
		Graphics3D[Flatten[{Hue[(#-First[tn])/(Last[tn]-First[tn])],lineopt,
		Line[Traj3D[Traj,t]/.t->#]}]
		],
		opt,AspectRatio->Automatic, DisplayFunction -> $DisplayFunction,
		PlotLabel->StringJoin["t=",ToString[N[#]]]]&,tn]
	];


(* ContourLines of a function in 3D *)

PlotContourLine3D[f_,{x_,x0_,x1_},{y_,y0_,y1_},{t_,t0_,t1_,dt_},opt___]:=
  Module[{a,u,s},
    Return[Show[
          Table[
            Graphics3D[
             First[Graphics[ContourPlot[Evaluate[f],{x,x0,x1},{y,y0,y1},opt,
                      ContourShading->False,
                      DisplayFunction->Identity]]
                  ]/.Line[u_]:>Line[Map[Prepend[#,t]&,u]]],{t,t0,t1,dt}],
        			DisplayFunction->$DisplayFunction,
				Axes->True]];]


(* ============================================================= *)

Other unsupported funtions

(* ============================================================= *)


(* PhasePlot Animation *)

AnimatePhasePlot[{xdot_,ydot_},Icond_,{t,t0_,T_},
		{pa_,pa0_,pa1_},params___]:=
	AnimatePhasePlot[{xdot,ydot},Icond,{t,t0,T},
				{pa,pa0,pa1,(pa1-pa0)/4},params];

AnimatePhasePlot[{xdot_,ydot_},Icond_,{t,t0_,T_},
				{pa_,pa0_,pa1_,pad_},params___]:=
Table[PhasePlot[
		Function[t,DeSol2DTab[xdot,ydot,
		Icond,t0,T]],t0,T,params,
		PlotLabel->"Step:"ToString[pa]],{pa,pa0,pa1,pad}];

(* PhaseMap and solutions in 3D space *)

PhaseVolSol[Traj_,t0_,t1_,Opt___]:=
	Show[SolPlot3D[Traj,t0,t1,
	DisplayFunction->Identity],
		Graphics3D[Line[Traj3D[Traj,t1]]],
		Graphics3D[Line[Traj3D[Traj,t0]]],
		PlotRange->All,
	DisplayFunction->$DisplayFunction,Opt];

(* Phasemaps in one plot *)

PhaseMapTab[Traj_,t0_,T_,dt_,opt___]:=
	Show[Graphics[Map[Line[Traj[#]]&,#]&[
					Table[t,{t,t0,T,dt}]]
				 ],
			PlotRange->All, AspectRatio->Automatic,
			Axes->True,
			PlotLabel->FontForm[ "The phasemaps",
			 {"Times-Bold", 12}],
			AxesLabel->{"x","y"},opt
		];

ListPhaseMapTab[Traj_,t0_,T_,dt_,opt___]:=
	ListPlot[
		Flatten[
			Map[Traj[#]&,#]&[
					Table[t,{t,t0,T,dt}]]
		,1],
			PlotRange->All, AspectRatio->Automatic,
			Axes->True,
			PlotLabel->FontForm[ "The phasemaps",
			 {"Times-Bold", 12}],
			AxesLabel->{"x","y"}
		,opt];

(* Phasemaps in graphic array *)

PhaseMapTab1[Traj_,t0_,T_,dt_,part_]:=
	Show[
		GraphicsArray[
		Partition[
			Map[Graphics[Line[Traj[#]],
					PlotRange->{{-1,1},{-1,1}},
					AspectRatio->Automatic,Axes->True,
					Ticks->None,
					PlotLabel->FontForm["t="ToString[N[#]],
					 {"Times-Bold", 10}]
				]&,#]&[Table[t,{t,t0,T,dt}]]
			,part]
			],
			PlotLabel->FontForm[ "The phasemaps", 
			{"Times-Bold", 12}]
		];



(* PhasePlots and PhaseMaps together *)

PhaseMapTr[Traj_,t0_,T_]:=
		Show[PhasePlot[Traj,t0,T],PhaseMap[Traj,T]];

(* Soln.-s and multiple phasemaps in 3D *)

PhaseVolTab[Traj_,t0_,T_,dt_,Opt___]:=Show[SolPlot3D[Traj,t0,T,
DisplayFunction->Identity],
Graphics3D[Map[Line[Traj3D[Traj,#]]&,#]&[Table[t,{t,t0,T,dt}]]],
		PlotRange->All, AspectRatio->Automatic,
		Axes->True,
		PlotLabel->FontForm[ "The solutions and phasemaps",
		{"Times-Bold", 12}],AxesLabel->{"t","x","y"},
		DisplayFunction->$DisplayFunction];
